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ABSTRACT 

Star formation in the inner few hundred pc of the Galactic bulge occurs in a flattened molecular layer called 
the central molecular zone (CMZ). Serabyn & Morris (1996) suggest that the star formation in the CMZ has been 
sustained for the lifetime of the Galaxy, and that the resulting agglomeration of stars formed in the CMZ has 
resulted in the prominent r~ 2 stellar density cusp at the Galactic center having about the same physical extent 
as the CMZ. This "central cusp" is somewhat less flat than the CMZ; thus the population of stars formed in the 
CMZ appears to have diffused out to larger latitudes. We hypothesize that such vertical diffusion is driven by the 
scattering of stars off the giant molecular clouds (GMC) in the CMZ, and perform numerical simulations of the 
scattering between stars and GMCs in the presence of the non-axisymmetric background potential. The simulation 
results show that the time scale for an initially flattened stellar population to achieve an aspect ratio of the observed 
OH/IR stars in the inner bulge, 1 to 2 Gyr, agrees well with the estimated age of those OH/IR stars. 

Subject headings: celestial mechanics, stellar dynamics — Galaxy: center — Galaxy: kinematics and dynamics 
— methods: numerical 



1. INTRODUCTION 



It has long been known that the central bulge of the Galaxy 
has a distinct stellar density profile which approximately fol- 
lows an r~ 2 law (Becklin & Neugebauer 1968). This r~ 2 cusp 
has often been assumed to be merely the innermost part of the 
elderly bulge (Evans & de Zeeuw 1994, among others). How- 
ever, the fact that the Central Molecular Zone (CMZ) has the 
same maximum identifiable extent in Galactic longitude as this 
stellar cusp (~ 200 pc) led Serabyn & Morris (1996) to pro- 
pose that star formation has been steadily occurring in the CMZ 
throughout the lifetime of the Galaxy. 

The CMZ is a massive reservoir of molecular clouds unparal- 
leled in the Galaxy. It is an active region of star formation with 
a number of very young stellar clusters and ample ionized gas 
occupying the central ~ 400pc x lOOpc region {I x b; see Fig- 
ure p. The clouds within the CMZ have relatively high density 
(n > 10 4 cm" 3 ), high volume filling factor (f > 0.1), and sig- 
nificantly elevated temperatures (30-200 K, typically ~ 70 K; 
Hiittemeister et al. 1993, among others), and their total mass is 
estimated at 5- 10 x 1O 7 M (Serabyn & Morris 1996), repre- 
senting roughly 10% of our Galaxy's neutral gas content. The 
inevitable angular momentum loss in an orbiting gas disk re- 
sults in gas inflow along the Galactic plane from the outer bulge 
to the CMZ. Serabyn & Morris enumerate the mechanisms for 
the angular momentum loss as follows: 1) shear viscosity in 
the differentially rotating gas disk, 2) shocks associated with a 
bar potential, 3) cloud-cloud collisions, 4) dynamical friction 
of giant molecular clouds by field stars (Stark et al. 1991), 5) 
magnetic field viscosity (Morris 1996), and 6) dilution of spe- 
cific angular momentum by stellar mass loss material raining 
down out of the slowly rotating Galactic bulge (Jenkins & Bin- 
neyl994). 

Both the stellar population (r~ 2 cusp) and the CMZ are con- 
siderably flattened along the Galactic plane. While the aspect 
ratio of the bulge is 1.6-1.7 (Kent, Dame, & Fazio 1991; Wei- 
land et al. 1994), those of the stellar population and the CMZ 
are - 2.2 (Catchpole, Whitelock, & Glass 1990) and ~ 4 (Ser- 
abyn & Morris 1996), respectively. The mass of the stellar 



population within the - 100-pc extent of the CMZ, ~ 10 9 M Q , 
falls near the lower end of the range predicted for the stellar 
mass emerging from the CMZ over the Galaxy's lifetime, and 
Serabyn & Morris consider this as additional support for the 
association between the stellar population and the CMZ. 

Here we focus our attention on the difference of the aspect 
ratios between the r~ 2 stellar population and the CMZ. In the 
context of the sustained star formation hypothesis, one may in- 
terpret the observation that the aspect ratio of the stellar popula- 
tion is smaller than the CMZ in terms of the vertical diffusion of 
stars formed in the flatter CMZ by gravitational perturbations. 

The vertical heating of the stars born in the Galactic disk 
due to gravitational perturbations has been studied by many au- 
thors, especially for heating by giant molecular clouds (GMCs). 
Spitzer & Schwarzschild (1951, 1953) first proposed that ran- 
dom encounters of the stars formed in a thin disk with mas- 
sive interstellar clouds would heat up the stellar populations 
(the presence of GMCs was not known until ~ 20 years later). 
They used a Fokker-Planck model to calculate the evolution of 
the velocity dispersions. Lacey (1984) performed similar cal- 
culations but with a consideration of vertical motions of stars. 
Villumsen (1983, 1985) approached the same problem by inte- 
grating the equation of motion for a tracer population of non- 
self-gravitating stars evolving in the gravitational fields of the 
fixed background (disk and halo) and of the GMCs. On the 
other hand, Barbanis & Woltjer (1967) proposed that transient 
spiral density waves could heat up the stellar populations. A 
series of papers (Carlberg & Sellwood 1985; Carlberg 1987; 
Sellwood & Kahn 1991) has explored this mechanism numer- 
ically. However, neither of the above two mechanisms alone 
was able to explain the observations. The models involving 
scattering by the GMCs predicted the ratio of stellar velocity 
dispersions perpendicular to the plane and towards the Galactic 
center, <j z /<jr, to be larger than that observed. They also pre- 
dicted that the velocity dispersion increases with time as ~ f 25 , 
which is substantially slower than implied by the observations, 
~ f 5 (Wielen 1977). Meanwhile, spiral structure was found to 
be highly inefficient in increasing a z . Thus, both mechanisms 
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FIG. 1. — Integrated intensity of 12 CO, J = 1-0 emission, from the AT&T Bell Labs survey of the inner Galactic bulge region (Uchida, Bally, & Morris, 
unpublished work). Contour lines were drawn at intervals of 250 K ■ kms~', starting from 250 K ■ kms~' . 



needed to be considered to account for the observations, and 
the study of the stellar heating under the combined influence of 
both GMCs and spiral structure was made by Jenkins & Binney 
(1990), whose simulations were found to agree with observa- 
tions of spiral structure material near the Sun. 

Contrary to the case of the Galactic disk, the kinematic evo- 
lution of a population of stars formed in the CMZ via heating 
processes has not been studied. The kinematic environments 
in the inner bulge differ from those in the disk mainly in three 
respects: 1) the background gravitational potential in the inner 
bulge is neither axisymmetric nor plane-parallel; 2) the pres- 
ence of density waves strong enough to significantly heat the 
stellar population seems unlikely (weakly damped modes may 
be present though; see, e.g., Weinberg 1994); and 3) the dynam- 
ical time scale in the inner bulge is considerably shorter than in 
the disk. Because of 1), we may not apply the analytic study 
of the kinematic evolution of the disk stars directly to the inner 
bulge. In the present paper, we hypothesize that the stellar pop- 
ulation formed in the CMZ diffuses mainly by scattering off the 
GMCs and we study the kinematic and spatial evolution of stars 
in a flattened configuration by following the trajectory of non- 
self-gravitating particles under the gravitational potential of the 
background (inner bulge and nucleus) and the GMCs, similarly 
to the work by Villumsen (1985) for the case of the Galactic 
disk. We choose to approach this problem numerically, in order 
to obtain more tangible results. 

If the scattering off the GMCs causes the stellar population 
to diffuse out from the inner bulge faster than the rate at which 
stars are formed in the CMZ, the resulting stellar population 
would have a significantly larger extent than the CMZ, instead 
of forming a configuration having a scale similar to the CMZ, 
such as the r~ 2 stellar population that we currently observe. On 
the other hand, if the diffusion is highly inefficient, the resulting 
stellar population from the sustained star formation in the CMZ 
would have a configuration very similar to that of the CMZ. 
Since the time scale of the stellar diffusion in the inner bulge is 
an important parameter, the calculation of the vertical evolution 
of the stellar population and its comparison with observations 
are essential for assessing our hypothesis. We compare our sim- 
ulation results with the distribution of OH/IR stars, one of few 



object categories in the inner bulge that currently provides ap- 
proximate information on the age. 

Our simulation models are described in § |2|. We show and 
discuss our results in §§ || and [|. Finally, our findings are sum- 
marized in § H 

2. MODELS 

As discussed in § 1, the background potential in the inner 
bulge is neither axisymmetric nor plane-parallel, which com- 
plicates any analytic approach for the kinematic evolution of 
stars within it. We therefore study the problem by directly inte- 
grating the orbital motion of non-self-gravitating test particles. 
The total integration time of our simulations is either 1 or 3 Gyr. 
Our model parameters discussed below are given in Table [jj 

2.1. GMCs 

Many groups have conducted radio observations of molecu- 
lar clouds in the inner bulge (see Morris 1997 for a list of such 
observations). Here, for the GMC parameters in our simula- 
tions, we adopt a CO line survey by Oka et al. (1998). They 
identified > 15 molecular clouds with a size of ~ 30 pc and 
obtained a total molecular mass in all clouds of ~ 4 x 10 7 M Q . 
We take 20 for the number of GMCs, N G mc, and 2.5 x 10 7 or 
5 x 10 7 M Q for the total mass in GMCs, M CM c- For the sake of 
simplicity, a Plummer model is adopted for the potential by the 
GMC: 

<W~- 2 o 1/2 , (1) 

where otgmc is the mass of a GMC, tcMC the core radius of the 
Plummer model, and d the distance to the GMC. Considering 
the size of the GMCs in the inner bulge, we use ccmc of 5 or 
10 pc. 

Since only the projected positional distribution and the 
line-of-sight velocity of the GMCs are available, the three- 
dimensional positional and velocity distribution of GMCs must 
be assumed. We assume that the R (galactocentric radius pro- 
jected onto the Galactic plane; note that we use r for the three- 
dimensional galactocentric radius and I for the projection of R 
onto the plane of the sky) distribution ranges from 10 to 150 pc 
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Table 1 

Model Parameters for Diffusion Simulations 
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Model 


(kms -1 kptT 1 ) 


(M ) 


Ot-GMC 


(pc) 


ct 


(pc 2 G-' 


Myr 2 ) 
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65 


2.5 x 10 7 




10 


-1.8 


3.1 > 


: 10 4 


0.8 
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65 


2.5 x 10 7 




10 


-1.8 


3.1 > 


! 10 4 


0.7 


3 


65 


2.5 x 10 7 




10 


-1.5 


4.5 > 


! 10 3 


0.8 


4 


65 


2.5 x 10 7 




10 


-1.8 


7.8 > 


! 10" 


0.8 


5 


50 


2.5 x 10 7 




10 


-1.8 


3.1 > 


! 10 4 


0.8 


6 


65 


5.0 x 10 7 




10 


-1.8 


3.1 > 


! 10 4 


0.8 


7 


65 


2.5 x 10 7 





10 


-1.8 


3.1 > 


! 10 4 


0.8 


8 


65 


2.5 x 10 7 


-1 


5 


-1.8 


3.1 > 


; 10 4 


0.8 



Note. — Model 1 is our standard model. Variations from model 1 are 
denoted by boldfaces. For all models, a VJcy = 20kms _1 , a v , z = 25kms _1 , 
Ngmc = 20, y = 0.9, and a = r = 1 pc. 



and that the surface number density varies as R aGMC . Binney 
et al. (1991) proposed that GMCs inside the 180-pc molec- 
ular ring move along the so-called X2 orbits (Contopoulos & 
Mertzanides 1977) and are gradually transported inward at a 
constant velocity due to angular momentum loss. Since the gas 
is thought to flow inwards along the Galactic plane, if the gas is 
transported inward while conserving its mass, oigmc is required 
to be -1. We here try cxgmc — ~ 1 an d 0. The inward migration 
takes place on a time scale much longer than that of the X2 or- 
bital motion, thus we neglect the inward migration and assume 
that the GMC motions in the plane follow X2 orbits. We also 
assume that the CMZ maintains the same properties over the 
entire computational time interval. 

The information on the motion of GMCs along the z-axis 
and its dependence on R are very limited. Therefore we as- 
sume that the vertical velocity dispersion of the GMCs, af MC , 
is constant over R and simply impose a sinusoidal vertical mo- 
tion to the GMCs with an amplitude Hqmc(R) that corresponds 
to the adopted cr™ c . Under the standard background potential 
adopted below, the X2 orbital motion with a™ c = 25kms _1 
gives Hqmc nearly proportional to R, with Hqmc — 25 pc at 
R = 150 pc. We use Hqmc = 25 (7?/150pc) pc in the simulations. 
For comparison, the 15 giant molecular clouds observed by Oka 
et al. (1998) have an average z-distance of 18 pc, and Oka et 
al. estimate a™ c to be > 32kms _1 . The initial positions of 
the GMCs are randomly chosen, but follow the assumed radial 
distribution, R aaMC , and GMCs are assumed not to interact with 
each other. In order to lessen the CPU burden, we pre-calculate 
the X2 orbit of each GMC and save it in a table. Then, when the 
equation of motion for the stars is integrated, the x-y position of 
each GMC at a given time is obtained from the table and the 
z position is obtained from the sinusoidal motion with period 
AH GMC /a™ c , and with off c = 25kms- 1 . 

2.2. Background Potential 
We model the potential of the inner bulge as 



$ = $ 1 + 



«o 



2+a 



x2+y - + Z - 

yt) 4 



1/2 



(2b) 



(r+r ) 



(2a) 



where r is the galactocentric radius and the softening parame- 
ter ro, which has negligible effect on the results, is adopted sim- 
ply to avoid numerical difficulties. The first term in the right- 
hand-side of equation ([la|) is a softened power-law ellipsoidal 
potential, and the second term, negligible at radii beyond a few 
parsecs, is to account for the contribution to the potential by the 
central black hole. One could obtain a more complicated po- 
tential model based on the the density model that matches the 
observed luminosity profile (e.g., one based on an ellipsoidal 
density distribution), but the simple potential above has been 
adopted here for the sake of efficient force calculation. 

The ellipsoidal potential is assumed to rotate with a pattern 
speed n p = 50 or 65kms~ 1 kpc~ 1 , following Menzies (1990), 
Binney et al. (1991, 1997), Zhao et al. (1996), and Dehnen 
(1999). Of course, the potential in the outer bulge will not 
be well described with the above formula, but few inner bulge 
stars will have enough kinetic energy to reach the outer bulge. 
The masses in the outer bulge and the Galactic disk presum- 
ably have a non-negligible contribution to the potential in the 
inner bulge, but such contributions may be understood to be al- 
ready included in parameters yo an d Zo in equation (pb[). While 
Becklin & Neugebauer (1968) found the IR luminosity dis- 
tribution for r < 25 pc to be ~ r" L8 , the kinematic study by 
Lindqvist, Habing, & Winnberg (1992a) implies a density dis- 
tribution in the 30pc < r < lOOpc region proportional to ~ r" 1 5 
(for lpc < r < lOOpc, the profile becomes ~ r" 20 ). We use 
a = -1.8 or -1.5. For both ciq and ro, we adopt 1 pc. The fit 
to the observations with the dust model of Spergel et al. (1996) 
gives a density aspect ratio of 1 .67 (Binney, Gerhard, & Spergel 
1997). We find that our potential model with zo = 0.8 gives a 
density aspect ratio of ~ 1.7. We thus choose zo = 0.8 and 0.7. 
For the set of parameters adopted here, the transition from X; 
orbits to X2 orbits, which Binney et al. (1991) propose to be 
responsible for the 180-pc molecular ring, forms at r ~ 180 pc 
only whenyo > 0.9. Thus we adopt yo = 0.9. Finally, we choose 
$0 such that $o(l +r/ao) 2+a corresponds to a mass distribution 
with a total mass inside r = 30 pc of 8 x 10 7 M© (all models ex- 
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cept model 4; Lindqvist et al. 1992a) or 2 x 1O 8 M (model 4; 
Genzel & Townes 1987). The equipotential and isodensity con- 
tours of the above model are shown in Figure [| The density 
distribution of a model with zo = 0.7 is slightly pinched at the 
z-axis (perpendicular to the Galactic plane). However, the cal- 
culation actually uses the derivative of the potential, and nei- 
ther the potential nor its derivative has such a feature. Thus, the 
pinch in the density will not significantly affect our results. 

2.3. Test Particles 

At every 10 pc between radii of 10 and 150 pc, 300 test par- 
ticles (stars) are initially distributed in the plane around a ring 
of that radius, making a total of 4500 stars per simulation. This 
discrete initial distribution allows us to apply one simulation 
result to initial power-law distributions of stars with various ex- 
ponents by differently weighting the stars initially placed at dif- 
ferent R. Simulation results are analyzed here for a* = -2 and 
— 1, where a* is the exponent of the initial, power-law surface 
density profile of stars in the plane (R a * is a surface density pro- 
file projected along the z-axis, not along the line-of-sight). The 
stars' initial azimuthal angles are random. Stars initially move 
on X2 orbits, and Gaussian random motions parallel and per- 
pendicular to the plane are imposed, with standard deviations 
o'v.xy = 10-20kms _1 and a ViZ = 12-25kms _1 , respectively. The 
value of 25kms _1 for a ViZ is chosen to match the a™ c value in- 
troduced earlier in deriving Hqmc = 25(/?/150pc) pc. Stars are 
first located at z = and then are propagated for approximately 
10 orbital periods to make the initial distribution phase-mixed. 
Since the dynamical relaxation time in the inner bulge is larger 
than the Hubble time, we neglect the interaction between stars. 

2.4. Numerical Method 

We have tested three numerical methods for integrating the 
equation of motion: the adaptive step size Runge-Kutta method, 
the Richardson extrapolation method (Bulirsch-Stoer method, 
Press et al. 1992), and the predictor-corrector method (vari- 
able time step leapfrog method; Hut, Makino, & McMillan 
1995), and found that the Bulirsch-Stoer method in a rotating 
frame is the most efficient for our problem. We also tried Stoer- 
mer's scheme, which is for second-order conservative equations 
like the equation of motion without a derivative on the right- 
hand-side (see, e.g., Press et al. 1992), in conjunction with the 
Bulirsch-Stoer method, but found it to be slightly less efficient 
than our choice above. In our simulations, the energy of a single 
star is conserved to an accuracy of better than 3 x 10~ 6 during 
the entire time interval followed. 

3. RESULTS 

3.1. Scale Heights and Velocity Dispersions 

Figure |] shows the 3 Gyr evolution of the dispersions of 
vertical stellar positions and velocities for model 1 . Stars ini- 
tially at 30 < R/pc < 70 and at 80 < R/pc < 120 are sepa- 
rately shown. Both vertical heights and velocities significantly 
increase in the first 1 Gyr. The initial a z values are smaller 
for the group of stars with smaller R because the adopted po- 
tential for the bulge has a larger absolute vertical gradient at 
smaller R. It is interesting that the increase of vertical velocity 
is almost independent of R, since stars closer to the Galactic 
center are expected to encounter GMCs more frequently than 
the ones at larger distances for this model (rotational velocity 
is nearly independent of R and the surface density of GMCs is 
larger at smaller R for this model) and thus one would expect 



the velocities of stars closer to the center to grow faster. We 
attribute this independence of vertical velocity evolution on R 
to the rapid evolution of the motions in the plane for stars at 
smaller R. Figure |] shows the average radii in the plane, R, for 
two radial groups of stars of model 1 . R of the stellar group with 
initially smaller R doubles and becomes comparable to that of 
the group with initially larger R in only 0.5 Gyr. Furthermore, 
we find that the average ratio of apocenter to pericenter radii 
of stellar orbits rapidly increases from ~ 1 .5 to ~ 4 in 0.5 Gyr 
for both groups, and keeps increasing afterwards. Once the or- 
bit of a star becomes highly eccentric, the star may encounter 
any GMC in the CMZ during its orbital motion, regardless of 
its initial location in the CMZ (or, regardless of its initial local 
number density of GMCs). 

Also shown in Figure Q is the evolution of the velocity dis- 
persion in the plane, cr v . xy , for model 1 . Unlike a v ~, er, uv ends 
its rapid increase and flattens after t ~ 1 Gyr. The radial veloc- 
ity dispersion in the plane, <j v ,r, and the tangential velocity in 
the plane, a V: g, have comparable values (0.9 < a Vt g /(T v ,r is 1-1)> 
implying that the velocity dispersion in the plane is nearly 
isotropic. We also find that cr VlZ /(T Vt R stays mostly between 0.7 
and 0.8, which is slightly larger than the value of 0.6 found by 
Villumsen (1985) for the Galactic disk case. 

The growth of the velocity dispersion may be described by 
the diffusion equation 



whose solution has a functional form of 

a(t) = a Q (l+'-y , (4) 

with n = 1/ (m + 2). When (>r, equation (^) simplifies to 

cr(f)ocf". (5) 

Villumsen (1985) found that the vertical velocity dispersion of 
stars in the Galactic disk is well fit by equation (M) with n = 0.38 
(m = 0.63) for the first 1 Gyr, and by equation (|j) with n = 0.31 
(m = 1.2) for t > 0.3 Gyr. These values of n are only slightly 
larger than the theoretical estimates by Lacey (1984), 0.33 and 
0.25, respectively. The theoretical estimation of these values 
for the inner Galactic bulge, however, is not easy because 1) 
the background potential is not plane-parallel, 2) the orbits of 
stars may not be approximated as epicyclic motions, and 3) the 
influence of GMCs may not be independent of each other be- 
cause of their large number density and softening radius, ccmc- 
Nonetheless, a comparison of the growth rate for the bulge and 
the disk may be instructive. We applied non-linear, 3-parameter 
least squares fits (eq. [Q]) to the early phase of our vertical 
velocity dispersion, and linear, 2-parameter least squares fits 
(eq. [|]|) to the later data, and find that the growth rates depend 
on the choice of time range to fit. Model 1 shows n = 0.30 
(m = 1 .33) for t < 1 Gyr and n = 0.20 (m = 3.0) for t > 1 Gyr, 
and model 3 gives n = 0.32 (m = 1.1) and n = 0.23 (m = 2.3), 
respectively. As in the Galactic disk, n becomes smaller as the 
scattering off the GMCs drives stars farther and farther away 
from where the GMCs reside. 

Next, we compare the evolution among different models. 
The velocity evolution of all our models is shown in Figure || 
Panel a) shows simulations with different bulge potential mod- 
els, b) with different bulge masses inside 30 pc and pattern 
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speeds, c) with different GMC masses, and d) with different 
GMC distributions and sizes. Panels a) and b) show that in 
general, the vertical velocity evolution is not sensitive to the pa- 
rameters used to model the bulge. The heating is only slightly 
more efficient when the potential is shallower (larger a). 

On the other hand, panels c) and d) of Figure || illustrate that 
the vertical velocity evolution is dependent on the mass and 
size of GMCs, but not on the distribution of GMCs. Stars are 
found to have highly eccentric orbits, which allows them to en- 
counter any GMC in the CMZ. Thus the heating rate is deter- 
mined mostly by the total number and mass of the GMCs in the 
CMZ, and not by the distribution of GMCs. 

Although not shown in Table [l] and Figure ||, we have tried 
other initial stellar velocity dispersions (<r w = lOkms -1 and 
<7 V , Z = 12kms _1 ), but found no significant difference from our 
standard model. The velocity dispersions rapidly increase by a 
factor of a few in the early phase, so unless the initial velocities 
are different by an order of magnitude, any initial difference in 
velocity dispersions that is smaller than the amount of increase 
in the early phase will soon be indistinguishable. 

3.2. Density Profiles and Aspect Ratios 

In this subsection, we discuss quantities with which we may 
more directly assess the hypotheses that the r~ 2 cluster is a re- 
sult of sustained star formation in the CMZ and that the GMCs 
are responsible for the vertical diffusion of newborn stars. To do 
so, we project the simulation stars onto the l-b plane of the sky 
and compare them with observations of a certain stellar type. 
From a variety of models of the inner Galaxy, the major axis 
of the bulge bar potential has been found to have an angle of 



15-30° to the east from the line connecting the Galactic center 
and the Sun (Binney et al. 1997, among others). Here we adopt 
20°. 

Our simulations initially have an equal number of stars for 
each linear R bin, making the exponent of the radial profile 
of the initial surface density of simulation stars in the plane, 
a*, equal to -1 (note that a* is different from a, the exponent 
for the bulge density profile). However, by properly weighting 
the stars in each bin, one may construct volume and projected 
(along the line-of-sight) surface density profiles (we denote the 
latter with £) for different a* values. Figure || shows the £ pro- 
file along the I and b axes at several epochs for a* = -2. It is in- 
teresting that £ along the / axis gradually decreases at I < 30 pc 
while maintaining its overall profile. The £ profile along the b 
axis experiences a significant flattening in its power-law slope. 

Plotted in Figure ffl is the evolution of the best-fit power-law 
exponent for the £ profile between 10 and 100 pc along the 
I and b axes. As anticipated from Figure ^, while the I axis 
slope does not vary much, that along the b axis continuously in- 
creases and approaches the former at the end of the simulation. 
The slope values for the a* = -2 case are steeper than those for 
the a* = -1 case by ~ 0.5 (when the distribution extends to in- 
finity, a difference of 1 in a* should result in a difference of 1 
in the power-law £ slopes; however, the stellar distribution in 
our simulations, as in reality, is finite, so the relation between 
a* and the £ slope is not as strong as in the infinite case). Us- 
ing a circular aperture, Becklin & Neugebauer (1968) obtained 
integrated near-infrared intensities of the inner Galactic bulge 
region as a function of aperture size and derived a surface in- 
tensity profile (projected along the line-of-sight) proportional to 
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FIG. 6. — Surface density profiles (projected along the line-of-sight) for model 1 along the / axis (a) and along the b axis (b) at four epochs, 0, 0.3, 1, and 3 Gyr. 
Stars that are less than 10 pc away from the axes are considered. The numbers in the plot denote the time in Gyr. S is in arbitrary units. The angle between the 
major axis of the bar and the Sun-Galactic center line is assumed to be 20°, and a* = -2 is adopted. The i-cr Poisson error in log S is smaller than 0.05 octave 
except for Log S < 0. 



r -o.8 enc i ose( j m ass in the 5 < r/pc < 100 region obtained 
by Lindqvist et al. (1992a) also implies a surface density pro- 
file of oc r"° 8 . On the other hand, by fitting an ellipse with an 
aspect ratio of 2.2 to the distribution of stars with dereddened, 
apparent K mag between 5 and 6 (Haller & Rieke 1999 esti- 
mate such stars to be a few 10 8 yr old), Catchpole et al. (1990) 
obtained a steeper exponent, —1.4, after correcting for the disk 
contamination, crowding, and dark clouds. The invariance of 
the /-axis slope over time in our simulations makes the /-axis 
slope indicative of its original slope. From Figure [7| we find 
that the deprojection of our simulation with a* = -1 approxi- 
mately reproduces the slope obtained by Becklin & Neugebauer 
and Lindqvist et al., and that with a» = -2 gives the slope ob- 
tained by Catchpole et al. In any case, it appears reasonable to 
assume that the time-averaged star formation efficiency in the 
inner bulge has a dependence between R~ 2 and R~ l in the plane. 

Finally, we compare the vertical thickening of our simulation 
stars with the observations. Since it is more instructive to com- 
pare with the observation of a coeval population (otherwise, the 
overlap of populations of different ages would blend out the 
age dependence of a certain quantity), we use the OH/IR stars 
observed by Lindqvist et al. (1992b) in the inner bulge (total 
number of 134). One of the measures of the thickness or flat- 
ness of a distribution is the aspect ratio of major (/) to minor 
(b) axes. However, the aspect ratio itself can be defined in sev- 
eral ways. Catchpole et al. (1990) measured the aspect ratio 
of the stellar distribution using observed isodensity contours. 
However, we find that a distribution of 1 34 stars is too sparse to 
obtain meaningful contours. For this reason, we determine the 
aspect ratio as the ratio of the major to minor axes of an ellipse 
which best fits the distribution of stars projected onto the plane 
of the sky. The best-fit ellipse was taken to be that for which 
the area per unit angle, taken as a function of angle on the sky 
measured from the center of the ellipse, was best matched by 
the function describing the number of stars per unit angle. The 
binning utilized to determine the latter function was chosen to 



ensure a statistically large number of stars in each angular bin. 
The orientation of the major axis of the fitted ellipse was fixed at 
the orientation of the Galactic plane. This is a one-dimensional 
fit, thus is less sensitive to the density profile, which is not a 
primary concern here. The observation by Lindqvist et al. was 
constrained to a cross-shaped area, but we limit the calculation 
of the aspect ratio of both observations and simulations to the 
stars in the 200 pc x 80pc area positioned at the Galactic center 
(the number of OH/IR stars observed by Lindqvist et al. in this 
area is 124). We here assume a» = -2, but we find that a* = -1 
gives nearly the same aspect ratios. 

Figure pi shows the evolution of the aspect ratios of all mod- 
els. The evolution is noticeably dependent only on the mass 
and size of GMCs as well as the bulge mass. The former de- 
pendence is a naturally expected phenomenon, and the latter 
dependence appears to be due to an initially larger aspect ratio 
of model 4, which has a larger vertical potential gradient, but 
the same initial vertical velocity dispersion as the other models. 

In spite of the dependence on the GMC parameters, all mod- 
els except model 4 (larger bulge mass case) evolve to having 
aspect ratio values below 3 in 1-2 Gyr. The same elliptical 
fit to the distribution of observed OH/IR stars in the central 
200pc x 80pc area gives an aspect ratio of 2. 6^' 5 (super- and 
subscripts denote 95 % confidence limits), shown in Figure |8| 
as a horizontal dashed line. Based on the observed bolometric 
magnitudes of the OH/IR stars (M ho i ~ -5; Jones et al. 1994 
& Blommaert et al. 1998) and the theoretical initial mass- 
bolometric magnitude relation of Vassiliadis & Wood (1993), 
Sjouwerman et al. (1999) estimates OH/IR stars in the inner 
bulge to be 1-2 Gyr old. Since the time scale for an initially 
flattened stellar population to achieve an aspect ratio of the ob- 
served OH/IR stars agrees well with the estimated age of the 
OH/IR stars, we conclude that scattering by the GMCs is one 
of the predominant mechanisms, perhaps the most important 
mechanism, for the vertical diffusion of stars in the inner Galac- 
tic bulge. 
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FIG. 7. — Evolution of the surface density slope (projected along the line-of-sightli'or model 1 along the / axis (solid lines) and along the b axis (dashed lines). 
The slopes are the best-fit power-law exponent to the distributions obtained in Figure nl The slopes are calculated for the a* = — 2 case (thick lines) and the a* = —1 
case (thin lines). 



It has been suggested by Baud et al. (1981) that the ex- 
pansion velocity of the circumstellar shell of an OH/IR star, 
Vexp, is a good statistical age indicator. For the Galactic inner 
bulge, Lindqvist et al. (1992a) confirmed this correlation by 
showing that the OH/IR stars there with v exp > 18kms _1 have 
younger kinematical and morphological characteristics than the 
ones with v exp < 18kms _1 . Our aspect ratio calculation for the 
same sample finds that such a division is clearer when the di- 
vision is made at 20kms _1 : while the aspect ratio of stars with 
Vexp < 20kms _1 is 2. S^\, that of stars with v exp > 20kms -1 is 
3.2^[-3. Although the latter, which is calculated from a sample 
of only 29, is rather uncertain, the rapid decrease of aspect ra- 
tios in our simulations supports the notion that the former group 
is younger than the latter group. 

Another measure for the comparison between our simula- 
tions and observations is the root-mean-square (rms) value of 
the Galactic height (Zrms)- This is also a projected quantity, but 
is sensitive to the vertical evolution of stars. We limit the cal- 
culation of Zrms to stars with \l\ < 20 pc and \b\ < 80 pc, in 
order to minimize the effects of planar evolution and of a few 
outliers with very large |z| values. The evolution of z rm s for 
all models is shown in Figure ^ along with the z rms value for 
the OH/IR stars in the same region, 22 pc (thick dashed lines). 
The figure is made assuming a* = -2, but we find that a* = -1 
typically gives only 10-20 % larger z„„ s values. All models, 
probably except model 4, reach the value for the OH/IR stars 
within ~ 2 Gyr, as in the case of the aspect ratio evolution. One 
noticeable behavior difference between Zrms an d the aspect ratio 
is its largely different convergent values between models 1 and 
6. Thus, Zrms could be useful in comparing simulations with the 
observed distribution of stars older than OH/IR stars. 

4. DISCUSSION 

The surface density profile along the I axis of model 1 at 
3 Gyr (Fig. |6^) shows a nearly flat slope near I = 30 pc. We 



find that this is caused by a relatively shallow volume density 
profile in the plane between r = 30 pc and 100 pc, p ~ r" 1 (see 
Fig. [Kj. In the case of model 3 (shallower bulge potential), 
the volume density profile becomes even shallower than r" 1 at 
3 Gyr. However, the profile of model 7, where the GMC distri- 
bution is uniform over R, does not change much over its simu- 
lation period, 1 Gyr. Thus the shallow density profiles shown 
in o>gmc = ~1 models are attributable to the relatively larger 
abundance of GMCs at smaller R, which more efficiently de- 
pletes stars in that region by heating. It is difficult to estimate 
the current volume density profiles in the inner bulge, both ob- 
servationally and from our simulations. Observationally, we 
only have projected surface density profiles and limited radial 
velocity information. To estimate the density profile from our 
simulations, one needs to assume a star formation history and a 
radial dependence of the star formation in the inner bulge. Until 
these are better understood, simulations may not give a strong 
constraint to the density profile. However, our simulations still 
seem to be able to suggest that the local density in the plane 
between 30 and 100 pc may be as shallow as or close to r" 1 
(compared to r" L8 obtained by Becklin & Neugebauer 1968 for 
r < 25 pc and r~ 2S by Lindqvist et al. 1992afor 1 < r/pc < 10). 
Such shallow profiles may carry important implications for star 
formation in the inner bulge, since a sphere (the inner bulge in 
this case) with a volume density profile close to r" 1 has negligi- 
ble tidal forces in it (the r" 1 density profile gives no tidal forces; 
if the density drop is shallower than r~ l , the tidal force even be- 
comes compressive). Thus if the local volume density in the 
inner bulge can be indeed shallower than, say, r" L5 , formation 
of stars there may not be as difficult as has been speculated (the 
galactic tidal field acts against the collapse or contraction of 
molecular clouds into a star). The enclosed mass as a function 
of galactocentric radius deduced by Lindqvist et al. (1992a; 
Fig. 10) shows a steep increase (thus shallow decrease in den- 
sity) between 30 and 100 pc (steeper than ~ r LS , implying a 
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FIG. 10. — Volume density profiles in the plane, p, for model 1 at four epochs, 0, 0.3, 1, and 3 Gyr. Stars that are less than 10 pc away from the plane are 
considered. The numbers in the plot denote the time in Gyr, and dashed lines indicate profiles of r _I and f 1 . p is in arbitrary units. a» = -2 is adopted. The \-a 
Poisson eiTor is smaller than 0.05 octave. 



density profile shallower than r ), and this supports our sim- 
ulation results that the density profile becomes shallower in the 
30 < r/pc < 100 region and also supports the above discussion 
on the possibly favorable star formation environment there. 

Considering the age of the Galaxy, OH/IR stars in the inner 
bulge are relatively young. Our comparison with a "young" 
population strongly supports our hypothesis that the GMCs are 
responsible for the vertical diffusion, but does not necessarily 
imply that the entire r~ 2 stellar population is a result of sus- 
tained star formation in the CMZ. To address this issue, one 
needs to account for the evolution and accumulation of stellar 
populations over the lifetime of the Galaxy. Since the evolution 
of aspect ratios in Figure [| appears to converge to an asymptotic 
value at times > 2 Gyr, we may expect that the agglomeration 
of stars formed in the CMZ over the lifetime of the Galaxy will 
have a stellar distribution not much different from that of the 
last stages of our simulations. A systematic observation that 
covers the full inner bulge region with good resolution and sen- 
sitivity is needed to take the next step in this analysis. This need 
may be fulfilled by the currently ongoing 2MASS survey. 

5. SUMMARY 

We have performed numerical simulations of scattering of 
young stars off the GMCs in the inner bulge. Several differ- 
ent bulge potentials, GMC models and distributions, and initial 
conditions of the stellar population have been considered. 

First we find that when clqmc = ~1> the stars inside ~ 70 pc 
experience rapid radial diffusion in the plane during the early 
phase. The consequences of this are that 1) later evolution of 
the inner stars is very similar to that of the outer stars (initial 
R > 80), and 2) the depletion of stars in the inner region makes 
the volume density profile in the plane between 30 and 100 pc 
quite shallow (possibly close to r~ l ). The latter is supported by 
the rapid increase in the enclosed mass between 30 and 100 pc 
(Lindqvist et al. 1992a), and implies that the tidal forces in 



that region may not be as hostile to star formation as has been 
previously conjectured. 

After a rapid initial rise during the first ~ 0.5 Gyr, especially 
for the inner stars, R linearly increases afterwards (in case of 
model 1, R reaches 150 pc in ~ 3 Gyr). The vertical positional 
dispersion of stars, on the other hand, grows almost linearly for 
the whole computational interval (model 1 reaches a z = 50 pc 
at - 3 Gyr). 

The projected surface density of test stars along the galactic 
plane inside 100 pc decreases as stars diffuse out both radially 
and vertically, but roughly maintains its initial slope over our 
simulation period. On the other hand, the slope of the projected 
surface density profile along the b (vertical) axis becomes sig- 
nificantly shallower during the same period and makes the over- 
all stellar population rounder. 

The comparison between observations and our simulations 
of the Z-axis surface density profile (projected along the line- 
of-sight) in the 10 < r/pc < 100 region suggests that newborn 
stars have an initial surface density profile (projected along the 
z-axis in the plane) of R~ 2 -R~ l . However, the estimation of 3- 
dimensional structure from 2-dimensional observation requires 
more careful study in the context of non-axisymmetric stellar 
models in order to reveal the true volume density profile in the 
inner bulge. 

The aspect ratios of our young stellar configuration are ini- 
tially larger than 6, but the ratios become smaller than 3 in 1- 
2 Gyr, except for one model. This is because the scattering 
off GMCs diffuses stars preferentially in the vertical direction 
(since the density gradient is greater in that direction), except 
in the very early phase. The aspect ratios converge to values 
between 1.5 and 2.5 in 1-2 Gyr. Since OH/IR stars in the inner 
bulge , which are estimated to be 1-2 Gyr old (Sjouwerman et 
al. 1999), have an aspect ratio of 2.6, we conclude that scatter- 
ing by the GMCs is one of the predominant mechanisms, and 
possibly the most important for the vertical diffusion of stars in 
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the inner Galactic bulge. 
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